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A statistical mechanical study of fiuidized granular media is presented. Using a special energy 
■ injection mechanism, homogeneous fiuidized stationary states are obtained. Molecular dynamics 

' simulations and theoretical analysis of the inelastic hard-disk model show that there is a large 

^SJ , asymmetry in the two-particle distribution function between pairs that approach and separate. 

Large velocity correlations appear in the postcoUisional states, due to the dissipative character of 
the collision rule. These correlations can be well characterized by a state dependent pair correlation 
, function at contact. It is also found that velocity correlations are present for pairs that are about 

' to collide. Particles arrive to collisions with a higher probability that their velocities are parallel 

' rather that antiparallel. These dynamical correlations lead to a decrease of the pressure and of the 

collision frequency as compared to their Enskog values. A phenomenological modified equation of 
state is presented. 
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I. INTRODUCTION 



\ There has been much interest recently devoted to the description of granular media. The understanding of their 
^ ■ physical properties is important because they appear in many different phenomena taking place in our daily life as well 
I ; as in various industrial processes. Experiments have been devised which have permitted to focus on many different 
. and new aspects. When subjected to injection of energy (vibrated plates for example), they present many similarities 
■ with fluids: convection take place, patterns can form and time-periodicity may be observed 
I . For those fiuidized granular media, careful investigations have permitted to test the adequacy of hydrodynamical 
' ^ continuum equations generalized to take into account the dissipation of energy due to the inelasticity of the collisions 
between the grains. Those equations are usually closed with the help of an equation of state and of phenomenological 
laws inspired by the hard-sphere model for fluids |9|-[l4[. Also, kinetic theory models have been studied obtaining 
good accord with computer simulations in the low density cases ■ 

In most of these works, granular media are described microscopically by means of the Inelastic Hard Sphere (IHS) 
' ^ ' model. Grains are modeled as soft hard spheres that dissipate energy at collision through a constant restitution 
, coefficient a. It is quite remarkable that such a simple model can already account for many peculiarities of the 
' observed behavior |18-20 



In this paper, we would like to go one step further: namely to investigate the IHS model as if it was a genuine 
statistical mechanical model for a granular fluid, and examine the effects of the inelasticity on the statistical properties 
of the fluid. One first obvious difference is the absence of any equilibrium state. Let to itself, an assembly of inelastic 
hard spheres will evolve towards a final state with no motion. One can however achieve stationary states by allowing 
contact with an energy reservoir, but those states are noncquilibrium states with a permanent energy flow coming 
from the reservoir and dissipated at collisions into the internal elasticity of the grains (and neglected for our purpose) . 
It is our aim to describe these Non Equilibrium Steady State (NESS), considering the dissipativity g, defined as 
q—{\ — a) 12^ as the parameter responsible for deviations from the equilibrium hard sphere system. 

Letting a system composed of IHS grains is let to evolve freely with periodic boundary conditions, it cools down 
homogeneously. This homogeneous coohng state (HCS) has been widely studied jl^jl^,^,^, and because of its 
O ' simplicity the HCS has been used as a reference state to build up theories for non-homogeneous states. It is the 
analog of the equilibrium state (with a Maxwellian distribution) for elastic systems p^-ps]]. Among other features, it 
has been shown that, in this state, the velocity distribution function is not Maxwellian, having a nonvanishing fourth 
cumulant and a long velocity tail ||2^ . Also, long range velocity correlations are developed, in the form of vortex 
^ , fluctuations p6| |. 

' It has been realized that, for any given value of the inelasticity, whatever small, there is always a size over which 
the homogeneous reference state looses its stability. The system undergoes a transition towards a state with a shear 
flow spontaneously developing and possibly towards an inhomogeneous clustered state |l^,|l^,|l^,|^|8) . In order to 
characterize intrinsic properties like an equation of state, one needs to restrict to consider stable homogeneous states 
before any instability occurs. This means in general to consider slight inelasticities and finite system sizes where 1/N 
effects are present, N being the number of grains. 

Computer simulation techniques has proven very valuable in granular studies: like in molecular dynamics (MD), 
computers are used to integrate numerically the equations of motion of a few hundred to a few thousand grains whose 
interaction is usually limited to inelastic collisions. The computed behavior compares generally well with experiments. 
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which give confidence in the adequacy of the inelastic hard sphere model for accounting for most of the specificities 
of granular fluids. We are going to use such a technique with a mechanism added in order to keep the kinetic energy 
of the grains constant. 

The essential result obtained is the characterization of short range correlations which are built in granular fluid 
models and which determine both the static (equation of state) and dynamical properties (kinetic equation) of the 
granular fluid. The pressure of the stationary state is computed from its mechanical definition and we can show that 
it is related to the pair correlation function at contact for precollisional configurations only. This induces an extra 
variable in the pair correlation function, the angle between relative distance and relative velocity at collision. 

Second, the kinetic description of the stationary state is well described by the (modified) hierarchy equations. In 
order to write a closed equation for the one-particle distribution function, like the Enskog equation for hard spheres, 
some assumptions are made on the two particle distribution function. In Enskog's equation, it is assumed that the 
particles that collide do not have velocity correlations (or dynamical correlations). When two particles are at contact 
with approaching velocities (precollisional state), it is assumed that the two particle distribution function /*^^^ can be 



where X = <7(c^) is the pair correlation function at contact and the Heaviside function guarantees that it is a prec- 
ollisional state. Note that in elastic systems at equilibrium, the previous relation is exact. We show here that this 
assumption is no more true for inelastic hard spheres, with an increased probability that colliding disks or spheres 
are parallel rather than antiparallel, at least for moderately dense fiuids. Such velocity correlations might also appear 
in elastic fluids but only under nonequilibrium regimes (typically over macroscopic scales). Here they appear on 
the shortest scales, those which correspond to molecular distances, only because of the dissipative character of the 
microscopic dynamics. 

The article is organized as follows: first, we present with some detail the simulation technique which has been 
used. As already mentioned, we simulate grain dynamics at constant energy: at every collision, all particle velocities 
are rescaled so as to keep the total kinetic energy constant. We present the technique and emphasize in particular 
its equivalence with a time-rescaling change. Next, we analyze the pressure and the pair-correlation function in the 
stationary state by making time averages of observables and configurations. The analysis lead in particular to the 
need of considering velocity correlations which, as will be shown, are also responsible for a pressure drop. This is 
done in the fourth section, before the conclusions. 



In the HCS the system is in a non-stationary state that does not allow to make time averages of different properties. 
It is then necessary to inject energy to the system in order to keep is stationary. Different methods have been used, 
mainly the injection of energy trough the walls (vibrating or stochastic), by external fields (like the flow in a pipe) or 
by stochastic forces acting on the particles (see for example ) . These methods have the disadvantage of destroying 
the homogeneity of the fluid or, in the case of the stochastic forces, adding other dynamical effects ]3l| ]. 

In this article we use a thermalizing method that both preserves the homogeneity and the dynamical properties of 
the granular fluid. Formally, each time two particles collide, the dissipated energy in reinjected immediately into the 
system by multiplying all particle velocities by the same factor. This factor is chosen each time in order to keep kinetic 
energy constant. As the IHS model does not have an intrinsic energy (or time) scale the rescaling of all velocities does 
not change the evolution of the system, the collision sequence (thus all physical phenomena) is preserved but viewed 
at a different speed. In this sense this thermostating mechanism is the most appropriate for the HCS. 

In this nonequilibrium steady state (NESS) all macroscopic quantities are stationary and can be averaged in 
numerical simulations. If we set the initial energy to give a temperature equal to one, then all averaged quantities 
correspond to computing them at this temperature. To obtain the value at another temperature, simple dimensional 
analysis gives the desired value. For example, if pness is the computed pressure in the NESS, the actual value of the 
pressure at another temperature T is p = pnessT. 

The IHS model is simulated using event driven molecular dynamics |32l . A direct implementation of the thermalizing 
method would lead to a computational cost proportional to the total number of particles for each collision, making 
impossible to make long simulations. But we can take advantage of the lack of an energy scale to make it much faster. 
In effect, rescaling all velocities is equivalent to rescaling the time. We deflne a rescaled time by the relation ds = jdt 
with 



well approximated by |Q 
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) = x/^'nvi)/('Hv2)e(-vi2 • ri2)<5(ri2 - a) 



(1) 



II. CONSTANT ENERGY SIMULATIONS OF THE IHS MODEL 
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where E is the total energy in the system. Note that 7 is a piecewise constant decreasing function. It can directly 
checked that if rescaled velocities are defined as v = dr/ds, then the rescaled kinetic energy {K — ^ ■y^'^) is conserved. 
This transformation have been successfully used to study the shearing instability ||2^ . 

Then, the simulation is done for the cooling IHS model (that is, no velocity rescaling is done) but at each collision 
the new kinetic energy is computed and 7 is evaluated using Eq. ^. In the simulation, then, the energy and the average 
velocity decrease, but the rescaled velocities are computed using v = v/7. Having the rescaled velocities, all the static 
properties can be computed directly. Since the function 7 is piecewise constant, the s-time can be integrated in the 
simulation; this allows to make periodical measurements in the NESS (equally spaced in s). 

Finally, to avoid roundoff errors, each time the kinetic energy has decreased by some orders of magnitude (typically 
10~^ of the initial value) a real velocity rescaling is done to put again the temperature equal to one. In this process 
also the center of mass velocity is subtracted because, otherwise, small errors will be amplified by the rescaling. 

As already mentioned the HCS becomes unstable (shearing instability) for certain values of dissipativity, density, 
and system size. Given values for the dissipativity and density, the system size (number of particles) is constrained 
to be smaller than a certain value to keep the system homogeneous. In 2D for a low density and low dissipation case, 
it is given by 

TT 

■^max — (3) 

qn 

For larger densities there are more complex expressions written in terms of the transport coefficients, but in all 
cases there is a maximum system size over which the system becomes unstable. This phenomenon does not allow to 
simulate large systems (for example, for n = 0.2 and q ~ 0.02, Nmax = 800) and finite size effects are obtained. To 
avoid these effects, simulations with different number of particles are done and then the results are extrapolated to 
the infinite size limit. The extrapolation is done using the standard model for any quantity that does not vanish in 
this limit, that is, for any quantity A its size dependence is modeled as 

A{N) = + Ai/N (4) 

When Njnax is smaller than around 1000, the effects of the shearing instability are present before the critical value 
as fluctuations in the k = 2Tr/L mode of the shearing velocity |^^. These fluctuations have both large amplitudes 
and long correlation times (divergent at the critical point) . Some quantities couple strongly with this fluctuating field 
(for example, the pressure) and small system sizes must be studied to extrapolate to infinity. 

We present series of simulations in 2D for the IHS model in the non equilibrium steady state already described. 
Dimensions are chosen such that the disk diameter, particle masses and initial temperature are set equal to one. The 
simulations are done at three different number densities {n = N/V): n = 0.05, n = 0.1, and n — 0.2. In each case the 
dissipation coefficient takes the values q — q — 0.002, q — 0.005, q — 0.01, and q — 0.02. Given the small values of 
the dissipation coefficients, long simulations are needed to obtain good statistics in order to identify the g-dependence 
of the different quantities. 

As the simulations are done for small dissipations, the results are presented as series in the coefficient q. Also, when 
possible, the results are condensed in power expansions in the density, but these expressions must not be interpreted 
as assuming that these quantities are necessarily analytic in density or dissipativity. 

Finally, unit s are chosen such that the disk diameter a and the particles masses are set to one. Also, the temperature 
in the NESS is fixed to one. 



III. VIRIAL PRESSURE 



As in granular media there is no equivalent to a free energy or an entropy, the pressure can not be defined thermo- 
dynamically but only mechanically. We use the virial expression for the pressure that is a mechanical definition 
valid for any isotropic system composed of particles that interact with pairwise forces. If the total volume is V and 
the total kinetic energy is K , in 2D the pressure is given by 




In the NESS, kinetic energy is constant and equal to N and the forces are impulsive ones. Using the collision rule, 
it is easy to check that the forces are given by 

Fy = - ( 1 - q) I Vy • h] I f,j 5{t~Uj) (6) 
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where is the collision time for the pair i—j, and the relative velocity is evaluated before the collision. 
Then, if the average is written as a time average, the delta functions are integrated giving 

colls 

where t is the averaging time, and it has been used that particles diameter is set to one. 
In simulations we measure the quantity 



colls 



colls 



in terms of which the total pressure is given by 

?7 = n + 

4 



nu , . 

p = n+ —pi (9) 



where v is the collision frequency. 

The results from the simulations, fitted linearly with q are 



pi{n = 


= 0.05) = 


1.772(1 - 


9) 


- 0.163g 


pi{n 


= 0.1) = 


1.772(1 - 


q) 


- 0.262g 


pi{n 


= 0.2) = 


1.772(1 - 


q) 


- 0.528g 



(10) 



That collected can be fitted to 

pi ^ ^{1 -q)- 2.67 qn (11) 
The average in Eq. || can also be written also as an ensemble average 

(l-'?).|;(2)(i,2)|vi2-ri2r 



p = n + 



W 

e(-vi2 •ri2)rfri(ivirfv2d6' (12) 

where 9 is the angle between the relative velocity and the relative position. Assuming lack of velocity correlations 

at collisions, the two particle distribution function can be replaced by Eq. |l|. In this case, the integrals can be done 
explicitly giving 

P = n+^^ML^ (13) 
Then, in this approximation, pi can be expressed as 

Pr = ^(l-.) (14) 



To study the validity of this approximation and compare with the numerical results (|l l|) , we need to study the pair 
correlation function at contact, X: a-nd the collision frequency, v. 



A. Pair correlation function at contact 



For hard particles systems, the pair correlation fmiction at contact x is defined in elastic systems as x — g{r=a^), 
where g(r) is the pair correlation function and a is the particle diameter. In granular media, this definition is somewhat 
ambiguous and a direct application of the classical computational methods |l3^ to obtain x does not give the most 
relevant result. 

The pair correlation function is defined in elastic systems as the probability of having two particles separated at 
a distance r, with an adequate normalization. This definition does not take into account the relative motion of the 
two particles since it is known that in equilibrium, positions are not correlated with velocities. But, as it is shown 
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bellow, in granular media positions and velocities are highly correlated even in the low density limit, having the pair 
correlation function different behavior if the two particles approach or separate. 

In Enskog's kinetic theory, the x factor that appears in Eq. in the coUisional term, in the virial pressure, and 
in the transport coefficients must be understood as the pair correlation function at contact for particles that are 
approaching and not for the ones that separate. In what follows we define and describe the basic properties of a 
generalized x coefficient that depends on the dynamical state of the particles. 

We define the pair correlation function g{r, 9) as proportional to the number of pairs that are separated by a distance 
r and there is an angle between the relative position and the relative velocity. The normalization is chosen such 
that g goes to one for large distances. This function is computed in an analogous way as the usual pair correlation 
function [Q. The generalized pair correlation function at contact is defined as x(^) = g{f = o"'',^). 

It can be identified two classes of pairs according to the value of the angle 9: if cos(9) is positive then the particles 
are separating (postcoUisional states), and if it is negative the two particles approach (precoUisional states). In the 
Appendix ^ we show that the postcoUisional part of x(6') can be expressed in terms of the precoUisional one . Using 
the collision rule and the conservation of probability during a collision it is found that if 9 < tt/2 (postcoUisional) 



X{9) = [cos(6l)2 + sm{9f] \[n^ tan-\atan{9))] (15) 

where [tt — tan~^(Q;tan(0))] > 7r/2 is the precoUisional angle that gives the postcoUisional angle 9. 

In the Enskog approximation, the precoUisional pair correlation function is a constant xo- In this approximation, 
the complete function is 

(ff.jxo cos(0)<O 
^ 1 Xo [cos(0)' + a'sin(0)2]-' cos(0) > ^ ^ 

This function is discontinuous at 9 = ±7r/2. Its average is xo^ (l + that can be understood knowing that the 
postcoUisional normal relative velocity is a times smaller than the precoUisional one. That makes that the colliding 
pair rests 1/a times longer in the postcoUisional position, giving the previously obtained average (formally the time 
the particles stay at contact is zero, but x is the limit of g{r) that can be well defined for finite bins where the pair 
stays finite times). 

This kind of discontinuity in x(^) has been found in sheared elastic fluids [^9| . In that case the origin of the 
discontinuity is the nonequilibrium character of the one-particle distribution function. In our case it is originated by 
the non-conservative collision rule. 

The discontinuity in x(^) makes that the extrapolation of g{r,9) to contact presents numerical problems. For the 
precoUisional angles g is a smooth function but for the postcoUisional ones there is a discontinuity line arriving to 
9 = tt/2, that prevents from obtaining good estimations of x for angles slightly below 7r/2. In Fig. Q it is presented 
a numerical estimation of x(^) obtained in MD simulations for a low density case. The comparison with the Enskog 
theoretical value (Eq. |l^) is good except for a 80° < 9 < 90° where is was expected to fail. 

In the clustering regime of large systems (where the g/^tem is no more homogeneous) it has been observed a large 
dependence of x on 9, for the precoUisional angles ||3^,^. In our case, statistical errors in the results of molecular 
dynamics simulations do not allow to determine if there is a dependence for precoUisional angles, contrary to the clear 
dependence for the postcoUisional angles. Nevertheless if x depends on the precoUisional angle, it is small and not 
like the one reported in the clustering regime. 

We define the pair correlation function at contact, Xi as the average of x(^) over the precoUisional angles only. The 
simulation results for x, fitted linearly with q, are 

X(n = 0.05) = 1 + 0.0655 + 0.051g (17) 
X{n = 0.1) = 1 + 0.1389 + 0.053g 
X(n = 0.2) = 1 + 0.3129 + 0.070g 

An empirical expression for x in the 2D elastic case is |3^ . 

7r(25 - imr) , , 

^i'l = 4(4-n.)> 

A comparison with the measured values show that, for the IHS model, there is a very small dependence on q. 
Within the statistical error it can be said that x does not depend on q and its value is the same as in the elastic case. 
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B. Velocity distribution and collision frequency 



It is known that the velocity distribution function for the HCS of the IHS model is not a Maxwellian but a distorted 
one. In the low velocity region, it is predicted that the fourth cumulant is not zero and using the Boltzmann-Enskog 
equation it has been predicted its value |^^. The fourth cumulant is defined as 

fc4 = ^ ' ' (19) 

and the predicted value using the Boltzmann-Enskog equation is for small values of q 

k^K -2q + 22.^1bq^ (20) 

The fourth cumulant is measured in the MD simulations obtaining values consistent with the theoretical predictions. 
The collision frequency v is defined as the average number of collisions per particle and unit s-time. So defined, v 
is a stationary quantity. The simulation results fitted linearly with q are 

v{n = 0.05) = 0.189 - 0.0042g (21) 
v{n ^ 0.1) = 0.403- 0.022q 
u{n = 0.2) = 0.929- 0.277q 

The collision frequency can be estimated using the the approximation that there are no velocity correlations at 
contact (Eq. 0) . Taking into account the distortion from the Maxwellian velocity distribution, the collision frequency 
is given by p2[ 

V = 2^nx (\ - ^) (22) 

where x must be understood as the precoUisional value. 

Both X s-nd the term (1 — ^4/32) have a positive dependence on 5, but the MD results ( pl| ) show a negative one. 
This discrepancy shows that the hypothesis of lack of velocity correlations is false and must modified. Also, when the 
numerical values for x (Eq jrzl) and v (Eq. |2^) are replaced in the approximation (|lj) for pi the predicted pressure is 
larger than the MD value (|1 Ij) . That is, using the approximation that the two particle distribution function can be 
factorized as in Eq. ^ the predicted pressure is larger than the observed one. 

IV. VELOCITY CORRELATIONS AT COLLISIONS 

Special MD measurements are done to study the source of the discrepancies in pressure and collision frequency. We 
compute numerically coUisional averages sensible to the presence of velocity correlations at contact for precoUisional 
states. 

For a system composed of particles that interact with hard core forces, the collision probability is given in 2D by, 

dPeoi/(l, 2) CX |V12 • &\ /(2)(1, 2)5(ri2 - (7)6 (-V12 • ri2) dG dvi dv2 rfn dV2 (23) 

where a is the particle diameter and is the Heaviside function that restrict the velocities to precoUisional states. 
CoUisional averages, defined as the average of any quantity at every collision in the system, can be computed using 
the previous probability distribution. 

<^)coll<-' =^y"^(l,2)/('Hl2)|vi2-^|e(-vi2-<T)dadvidv2 (24) 

= — /^(l,2)/(2)(12)|v2-vi|dfedvidv2 (25) 
nv J 

where b is the impact parameter. The sign (— ) in the average means that A is evaluated with the precoUisional 
velocities and that r2 ri — a. It must be remarked that the above defined coUisional average only takes into account 
the precoUisional part of /^^^ 
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In the hypothesis of absence of precolhsional velocity correlations at contact, the two particle distribution function 
is written as in Eq. |l|. As the mean velocity vanishes, this approximation implies that the following coUisional average 
should vanish. 

(26) 



|V2 - Vi|/coll(-) 

It must be remarked that in an elastic fluid at equilibrium Eq. ^ is exact and F must vanish exactly. 

This property makes F a good test for the hypothesis of the lack of velocity correlations in stationary states. This 
quantity couples strongly with the velocity fluctuation so careful extrapolation to the infinite system must be done. 
In Fig. we show the extrapolation procedure for a typical series of simulations. For small system sizes F is negative 
as a consequence of the conservation of total momentum. In finite systems if one particle has a velocity vq the others 
have in average a velocity equal to —vq/{N — 1), leading to negative values for F. For large systems, the shearing 
instability appears giving rise to a dramatic increase of F. To get the extrapolated value to infinite system size, we 
consider systems up to a certain size where no signal of the instability is present {N = 1500 in the case of the figure). 

After extrapolation to infinity system size, the obtained results fitted linearly with q are 

(27) 



r{n = 


= 0.05) 


= 0.097(7 


r{n 


= 0.1) 


= 0.269(7 


F(n 


= 0.2) 


= 0.442(7 



That can be collected in the general expression 

T = 2.29 qn (28) 

This result means that in effect there are short range velocity correlations with origin in the dissipative character of 
the fluid. The fact that the correlation is positive means that particles arrive at collisions with velocities more parallel 
than if there were no correlations. This phenomenon can be understood in terms of recoUisions since, for the IHS 
model, the velocities of the particles after a collision become more parallel than in the elastic case. This parallelization 
has two effects: first it increases the probability of having a recollision (that is, after colliding with a third particle the 
tagged pair recoUides) and, second, when the pair recoUides, their velocities are correlated, being more parallel. The q 
dependence and density dependence can be well understood in this model because the parallelization is proportional 
to q and the recollision probability to n. 

This interpretation is also supported by the results in pressure and collision frequency. In effect, the parallelization 
of velocities after collisions produces that, in recoUisions, the approaching relative velocities are smaller, decreasing 
the collision frequency. Also the transfered momentum -which, when averaged, gives pi- is smaller at each collision. 
The combination of these two effects gives that the velocity correlations produce smaller values for the pressure. The 
fact that the effects of the correlations in F and pi are numerically similar is also an indication that the two effects 
have the same origin. 

At low density, precoUision velocity correlations are small. This explains the good agreement between the simulations 
and the theoretical predictions using Enskog's theory for x(^) in Fig- 

Precolhsional velocity correlations are also present in elastic fluids, but only under nonequilibrium conditions. Low 
frequency and long wavelength phenomena (where hydrodynamics is valid) are usually studied as perturbations over 
local equilibrium states where no velocity correlations exist. In granular fluids, the system is always out of equilibrium 
and velocity correlations are always present. For a given dissipation, there is no regime with no velocity correlations 
around which perturbation analysis can be done. 

In elastic fluids, the factorization (^ is exact at equilibrium and it allows to compute static properties like the 
pressure and the collision frequency. From this starting point, Enskog equation is built to describe the time evolution 
of the system in an approximate way (it neglects precolhsional velocity correlations in every regime) . To mimic the 
Enskog approach for granular fluids, we should take an equation that includes velocity correlations, even in the HCS, 
in order to predict accurately the static properties: pressure, collision frequency, dissipation rate. 

At low density and/or dissipation, the velocity correlations are small. Then, the Boltzmann-Enskog equation can 
be used to describe granular fluids with the same degree of approximation as for elastic fluids. For dense dissipative 
systems more complex theories, that take into account recoUisions (for example, Ring Equations |^,|l^), are necessary. 

V. CONCLUSIONS 

Different properties for the IHS model for granular fluids, put in the homogeneous cooling state, have been carefully 
studied. Using a time rescaling formalism it was possible to obtain precise averages in MD simulations, allowing to 
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study the dissipation and density dependence of these properties. Dissipation took values from q = to q = 0.02 and 
number density from n = 0.05 to n = 0.2. 

In all cases it was found that the velocity distribution is not Maxwellian and the fourth cumulant is different 
from zero. Its value does not depend on density and it is in good accord with the kinetic theory predictions. 

A distinction is made between the pair correlation function at contact for prccoUisional and postcoUisional states. 
The first is the one used in kinetic theory (Enskog's theory) and it was found that is has a very small dependence 
on dissipation, taking the same value than for elastic disks. The postcoUisional pair correlation function takes larger 
values and it can be fully predicted in terms of the precollisional function. 

Collisional averages indicate that particles that are about to collide are correlated in a non trivial way, particles 
arrive to collisions with velocities that are more parallel that in a elastic fluid. The computed correlation, that is 
proportional to density and dissipation, has its origin in recollisions: due to dissipation, particles that collide emerge 
with more parallel velocities than in the elastic case and, when they recoUide, their velocities are still more parallel. 
Results obtained for pressure and collision frequency also show the signature of velocity correlations at collisions. 
The effect of these correlations is to reduce the collision frequency and the transfered momentum at collisions, thus 
reducing the virial pressure. 

In elastic systems, velocity correlations are also present, but only in nonequilibrium regimes. The intensity of the 
correlations reduces as the system approaches equilibrium. In granular fluids, on the other hand, the dissipative char- 
acter of collisions puts the system always out of equilibrium, creating velocity correlations. The observed correlations 
are intrinsic to granular fluids since they are present in every regine. There is no need for special initial conditions 
or boundary conditions to obtain and compute them. This allowed us to compute them in a very simple regime, the 
homogeneous cooling state, with very high precision at the shortest possible scale, the microscopic one. 

The presence of these correlations implies that the Enskog factorization is insufficient to compute the static 
properties of the fluid: pressure, collision frequency, and dissipation rate. For elastic fluids, Enskog's equation, even 
if it is an approximation, accurately predicts static properties. An equivalent approach for dissipative systems would 
need the use of a kinetic theory that includes velocity correlations, even in the HCS. More complex theories like 
ring kinetic theory or mode coupling theories |39| are then needed to describe dense granular fluids at finite 

dissipation. 
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APPENDIX A: RELATION BETWEEN THE POST AND PRECOLLISIONAL PART OF x{9) 

In this appendix we will deduce the expression (|l^) . The deduction is based on the transformation of the distribution 
function at collisions and in geometrical aspects of the collision rule. 

The instantaneous character of binary collisions in the hard sphere (disk) system, implies that the two particle 
distribution function can be written as |2^ 

2)5(ri2 - a) - e(-ri2 • vi2)/f (1, 2) + ^e(ri2 • vi2)6*/,i"(l, 2) (Al) 

(2) 

The first term is the precollisional distribution function /q . The second term represent the postcoUisional distri- 
bution function, written in terms of the precollisional one. The operator h* has the effect of replacing the velocities 
with the precollisional values, and the factor comes from the the change in relative velocity and the Jacobian of 
the transformation [ p7[ |. 

The postcoUisional part of the pair distribution function is then 

/(2)(l,2)e(ri2 • vi2)5(ri2 - a) - a" Vo'^ (1*, 2*)e(ri2 • Vi2),5(ri2 - a) (A2) 

where 1* and 2* represent the state of the particles with precoUision velocities. 

To simplify notation we will consider a collision in 2D; for the 3D case, the analysis is similar and the results are 
summarized at the end. The geometry of the collision is represented in Fig. The postcoUisional relative velocity 
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is Vi2, the precoUisional relative velocity is vl2, and a is the vector that joins the centers of the two particles. The 
angles di (precoUision) and 6*2 (postcollision) are defined as 



'12 



The collision rule implies that 



^12 

a -1 /vi2 • cr 

O2 = COS 

tan(6'2) = -a"Uan(6'i) 



thus 



61 = TT — tan ^ (a tan(6'2)) 
The pair correlation function at contact x(^) is defined as 

-1 / V12 • cr 



cos 



Vl2 



dvi dv2 da 



(A3) 
(A4) 

(A5) 
(A6) 

(A7) 



where the factor guarantees the correct normalization. 

For postcoUisional angles, /'•^■'(l, 2) can be expressed in terms of the precoUisional velocities using (A2). Changing 
integration variables we obtain 



/f (r,2*) 6 



^1 / V12 ■ q- 

V ^^12 



dv\ ctv*2 da 



< tt/2 



(A8) 



where it has been used that dvi dv2 — adv^ Wl 



The argument of the delta function can be changed to precoUisional velocities using (A3), ( A4),and ( A6) 



x{e)^ — [cos\e)+a^sixv\9)\ 



— tan "'^(atan(0)) — cos ^ 



'12 ' • 



dv^ dvj da 6 < Tr/2 



(A9) 



where the transformation rule for delta function has been used. The integral can be identified as the pair correlation 
function at contact for the precoUisional angle [tt — tan^^ (a tan((?))] 



X{0) ^ [cos^ie)+a^sm^{e)] \ [n ~ tan-^ {a t£in{e))] e<7r/2 



(AlO) 



that is, the postcoUisional part of x(^) can be computed using the precolisional values. 

In 3D, we define the pair correlation function that depends on the solid angle Q that form the relative velocity and 
the relative position, where Ct is represented as usual by the angles and (f>. As the tangential components of the 
relative velocity are pres erved at the collision, 0i = 4>2- The change on the normal component of the relative velocity 
implies the relation ([A^). Using the generic relation of the delta function (5(57 — tl') — 5(6 — 9')5{(j) — svii{9)\ and 
that I sin(0i)| = | sin(6'2)|, it is found by a similar analysis as in the 2D case that 



X{n) = [cos^{e) + c? sin2(6')] ^ x{^*) d < tt/2 
where Ct* is the precoUisional solid angle. 



(All) 
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FIG. 1. Estimation of x(^) obtained in MD simulations (dots) compared to the theoretical prediction. The simulation 
parameters are N — 1000, n = 0.02, and q = 0.1. The discrepancy near 9 = it/2 is explained in the text. 




FIG. 2. Values of T obtained in MD simulations for an elastic case and a dissipative case as a function of the inverse of 
the number of particles N. The solid circle indicates the extrapolated value for the dissipative system. The global density is 
n = 0.1. 
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FIG. 3. Geometry of in inelastic collision. Vi2 (V12) is the incoming (outgoing) relative velocity and o is the normal vector 
to the collision. 
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